

library(tidyverse)
library(fixest)
library(broom)
library(ggthemes)
library(interflex)

#### Split samples for heterogeneous effects ####

# Democracy
models_democracy_annual <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                            scaled_elec_renew, scaled_fit, scaled_quota,
                            scaled_delegates, scaled_gcg,
                            scaled_cf_total_provided, scaled_cf_principal_provided,
                            scaled_cf_total_received, scaled_cf_principal_received)
                          ~ annual_average_temp_lag1 + annual_W_lag1 |
                            iso3c[year] + year, # trends
                          data = filter(analysis_df, democracy_dummy == 1)) %>% 
  as.list()
models_democracy_season <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                                   scaled_elec_renew, scaled_fit, scaled_quota,
                                   scaled_delegates, scaled_gcg,
                                   scaled_cf_total_provided, scaled_cf_principal_provided,
                                   scaled_cf_total_received, scaled_cf_principal_received)
                                 ~ hottest_season_temp_lag1 + season_W_lag1 |
                                   iso3c[year] + year, # trends
                                 data = filter(analysis_df, democracy_dummy == 1)) %>% 
  as.list()

# Autocracy
models_autocracy_annual <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                                   scaled_elec_renew, scaled_fit, scaled_quota,
                                   scaled_delegates, scaled_gcg,
                                   scaled_cf_total_provided, scaled_cf_principal_provided,
                                   scaled_cf_total_received, scaled_cf_principal_received)
                                 ~ annual_average_temp_lag1 + annual_W_lag1 |
                                   iso3c[year] + year, # trends
                                 data = filter(analysis_df, democracy_dummy == 0)) %>% 
  as.list()
models_autocracy_season <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                                   scaled_elec_renew, scaled_fit, scaled_quota,
                                   scaled_delegates, scaled_gcg,
                                   scaled_cf_total_provided, scaled_cf_principal_provided,
                                   scaled_cf_total_received, scaled_cf_principal_received)
                                 ~ hottest_season_temp_lag1 + season_W_lag1 |
                                   iso3c[year] + year, # trends
                                 data = filter(analysis_df, democracy_dummy == 0)) %>% 
  as.list()

# Rich
models_rich_annual <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                       scaled_elec_renew, scaled_fit, scaled_quota,
                       scaled_delegates, scaled_gcg,
                       scaled_cf_total_provided, scaled_cf_principal_provided,
                       scaled_cf_total_received, scaled_cf_principal_received)
                     ~ annual_average_temp_lag1 + annual_W_lag1 |
                       iso3c[year] + year, # trends
                     data = filter(analysis_df, rich_dummy == 1)) %>% 
  as.list()
models_rich_season <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                       scaled_elec_renew, scaled_fit, scaled_quota,
                       scaled_delegates, scaled_gcg,
                       scaled_cf_total_provided, scaled_cf_principal_provided,
                       scaled_cf_total_received, scaled_cf_principal_received)
                     ~ hottest_season_temp_lag1 + season_W_lag1 |
                       iso3c[year] + year, # trends
                     data = filter(analysis_df, rich_dummy == 1)) %>% 
  as.list()

# Poor 
models_poor_annual <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                       scaled_elec_renew, scaled_fit, scaled_quota,
                       scaled_delegates, scaled_gcg,
                       scaled_cf_total_provided, scaled_cf_principal_provided,
                       scaled_cf_total_received, scaled_cf_principal_received)
                       ~ annual_average_temp_lag1 + annual_W_lag1 |
                       iso3c[year] + year, # trends
                     data = filter(analysis_df, rich_dummy == 0)) %>% 
  as.list()
models_poor_season <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                              scaled_elec_renew, scaled_fit, scaled_quota,
                              scaled_delegates, scaled_gcg,
                              scaled_cf_total_provided, scaled_cf_principal_provided,
                              scaled_cf_total_received, scaled_cf_principal_received)
                            ~ hottest_season_temp_lag1 + season_W_lag1 |
                              iso3c[year] + year, # trends
                            data = filter(analysis_df, rich_dummy == 0)) %>% 
  as.list()

hte_list <- do.call(c, list(models_democracy_annual, models_democracy_season,
                            models_autocracy_annual, models_autocracy_season,
                            models_rich_annual, models_rich_season,
                            models_poor_annual, models_poor_season))

hte_results <-  setNames(as.data.frame(matrix(NA, nrow = 88, ncol = 10)), 
                               c("outcome", "treatment", "beta", "std_error",
                                 "t_stat", "p_value", "nobs", "n_countries",
                                 "min_year", "max_year"))

for (i in 1:88){
  hte_results[i, 1] <- hte_list[[i]]$fml %>% gsub("\\~", "", .) %>% .[2] # outcome
  hte_results[i, 2] <- hte_list[[i]]$fml %>% gsub("\\~", "", .) %>% .[3] # treatment
  hte_results[i, 3] <- hte_list[[i]] %>% tidy() %>% filter(!grepl("_W_",term)) %>% select(2) # beta
  hte_results[i, 4] <- hte_list[[i]] %>% tidy() %>% filter(!grepl("_W_",term)) %>% select(3) # std error
  hte_results[i, 5] <- hte_list[[i]] %>% tidy() %>% filter(!grepl("_W_",term)) %>% select(4) # t-stat
  hte_results[i, 6] <- hte_list[[i]] %>% tidy() %>% filter(!grepl("_W_",term)) %>% select(5) # p-val
  hte_results[i, 7] <- hte_list[[i]]$nobs 
  hte_results[i, 8] <- length(unique(hte_list[[i]]$fixef_id$iso3c))
  hte_results[i, 9] <- min(attr(hte_list[[i]]$fixef_id$year, "fixef_names"))
  hte_results[i, 10] <- max(attr(hte_list[[i]]$fixef_id$year, "fixef_names"))
}
hte_results$stratifier <- c(rep("Democracy", times = 11*2),
                            rep("Autocracy", times = 11*2),
                            rep("Rich", times = 11*2),
                            rep("Poor", times = 11*2))
hte_results$treatment <- str_sub(hte_results$treatment, start = 1L, end = -17L)



labels_vector <- hte_results %>% 
  filter(stratifier=="Democracy" & treatment == "annual_average_temp_lag1") %>% 
  mutate(labels = case_when(outcome=="scaled_climate_laws" ~ "Climate laws",
                            outcome=="scaled_fit" ~ "Feed-in tariff",
                            outcome=="scaled_quota" ~ "Renewables quota",
                            outcome=="scaled_elec_renew" ~ "Renewable electricity\n generation",
                            outcome=="scaled_ecp_zeros" ~ "Emissions-weighted\n carbon price",
                            outcome=="scaled_delegates" ~ "COP delegates",
                            outcome=="scaled_gcg" ~ "Climate institutional\n memberships",
                            outcome=="scaled_cf_total_provided" ~ "Climate finance\n provided",
                            outcome=="scaled_cf_principal_provided" ~ "Climate finance\n (principal) provided",
                            outcome=="scaled_cf_total_received" ~ "Climate finance\n received",
                            outcome=="scaled_cf_principal_received" ~ "Climate finance\n (principal) received")) %>%
  select(labels)


# Plot in four facets
rbind(hte_results) %>% 
  # Remove estimates from samples with fewer than 5 countries
  mutate(beta = ifelse(n_countries < 5, NA, beta),
         std_error = ifelse(n_countries < 5, NA, std_error)) %>% 
  group_by(treatment, stratifier) %>% 
  mutate(n = row_number()) %>% 
  ungroup() %>% 
  group_by(outcome, stratifier) %>% 
  mutate(m = row_number()) %>% 
  ungroup() %>% 
  mutate(index = ifelse(m == 1, n - 0.15, n + 0.15)) %>% 
  # Create point ID
  mutate(Temperature = case_when(str_detect(treatment, 'season')==TRUE ~ "Hottest season",
                                 str_detect(treatment, 'annual')==TRUE ~ "Annual")) %>% 
  # Create facet ID
  mutate(Facet = stratifier) %>% 
  # Get confidence intervals
  mutate(ci_lower = beta - 1.96*std_error,
         ci_lower_correction = beta - 2.84*std_error,
         ci_upper = beta + 1.96*std_error,
         ci_upper_correction = beta + 2.84*std_error) %>% 
  ggplot(., aes(beta, index, shape = Temperature)) +
  facet_wrap(~ Facet, nrow = 2) +
  coord_cartesian(xlim = c(-0.35,0.35)) +
  scale_x_continuous(breaks = c(-0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3)) +
  geom_segment(aes(x = ci_lower_correction, xend = ci_upper_correction,
                   y = index, yend = index),
               size = 1, color = "#79b321") +
  geom_segment(aes(x = ci_lower, xend = ci_upper, 
                   y = index, yend = index), 
               size = 2, color = "#456613") +
  # geom_point(aes(shape = Temperature), size = 2, color = "black") + 
  geom_point(size = 3.5, aes(shape = Temperature)) +
  scale_shape_manual(breaks=c("Annual","Hottest season"),
                     values = c(16,3)) +
  scale_color_manual(values = c("black","black")) +
  geom_vline(xintercept = 0, linetype = 1, colour = "black") +
  scale_y_continuous(trans = "reverse",
                     breaks = c(1:11),
                     labels = unique(labels_vector$labels)) +
  labs(x="Regression coefficient and confidence intervals", 
       y="") +
  geom_vline(xintercept = 0) +
  theme(panel.spacing = unit(1, "lines")) +
  theme_tufte(base_family = "Helvetica") +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill=NA),
        panel.grid.major.x = element_line(colour = "grey90"),
        panel.grid.minor.x = element_line(colour = "grey90"),
        # panel.grid.major.y = element_line(colour = "grey90"),
        panel.ontop = FALSE,
        text=element_text(size=16))
ggsave("text/coefplot_heterogeneous_trends.pdf", units = "in", height= 10, width = 10)


#### Continuous measures of democracy ####

# interflex(Y = "scaled_gcg",
#           D = "hottest_season_temp_lag1", 
#           X = "gdp_pcmkt", 
#           FE = c("iso3c", "year"),
#           estimator = "binning",
#           wald = F,
#           vcov.type = "cluster",
#           cl = "iso3c",
#           data = analysis_df, na.rm=T)

## Annual average temperature * GDP per capita (logged)
ga1 <- interflex(Y = "scaled_climate_laws",
                D = "annual_average_temp_lag1", 
                X = "gdp_pcmkt", 
                Z = "annual_W_lag1",
                FE = c("iso3c", "year"),
                estimator = "binning",
                wald = F,
                vcov.type = "cluster",
                cl = "iso3c",
                data = analysis_df, na.rm=T)
ga2 <- interflex(Y = "scaled_ecp_zeros",
                 D = "annual_average_temp_lag1", 
                 X = "gdp_pcmkt", 
                 Z = "annual_W_lag1",
                 FE = c("iso3c", "year"),
                estimator = "binning",
                wald = F,
                vcov.type = "cluster",
                cl = "iso3c",
                data = analysis_df, na.rm=T)
ga3 <- interflex(Y = "scaled_elec_renew",
                 D = "annual_average_temp_lag1", 
                 X = "gdp_pcmkt", 
                 Z = "annual_W_lag1",
                 FE = c("iso3c", "year"),
                estimator = "binning",
                wald = F,
                vcov.type = "cluster",
                cl = "iso3c",
                data = analysis_df, na.rm=T)
ga4 <- interflex(Y = "fit",
                 D = "annual_average_temp_lag1", 
                 X = "gdp_pcmkt", 
                 Z = "annual_W_lag1",
                 FE = c("iso3c", "year"),
                estimator = "binning",
                wald = F,
                vcov.type = "cluster",
                cl = "iso3c",
                data = analysis_df, na.rm=T)
ga5 <- interflex(Y = "quota",
                 D = "annual_average_temp_lag1", 
                 X = "gdp_pcmkt", 
                 Z = "annual_W_lag1",
                 FE = c("iso3c", "year"),
                estimator = "binning",
                wald = F,
                vcov.type = "cluster",
                cl = "iso3c",
                data = analysis_df, na.rm=T)
ga6 <- interflex(Y = "scaled_delegates",
                 D = "annual_average_temp_lag1", 
                 X = "gdp_pcmkt", 
                 Z = "annual_W_lag1",
                 FE = c("iso3c", "year"),
                estimator = "binning",
                wald = F,
                vcov.type = "cluster",
                cl = "iso3c",
                data = analysis_df, na.rm=T)
ga7 <- interflex(Y = "scaled_gcg",
                 D = "annual_average_temp_lag1", 
                 X = "gdp_pcmkt", 
                 Z = "annual_W_lag1",
                 FE = c("iso3c", "year"),
                estimator = "binning",
                wald = F,
                vcov.type = "cluster",
                cl = "iso3c",
                data = analysis_df, na.rm=T)
ga9 <- interflex(Y = "scaled_cf_principal_provided",
                 D = "annual_average_temp_lag1", 
                 X = "gdp_pcmkt", 
                 Z = "annual_W_lag1",
                 FE = c("iso3c", "year"),
                estimator = "binning",
                wald = F,
                vcov.type = "cluster",
                cl = "iso3c",
                data = analysis_df, na.rm=T)
ga8 <- interflex(Y = "scaled_cf_total_provided",
                 D = "annual_average_temp_lag1", 
                 X = "gdp_pcmkt", 
                 Z = "annual_W_lag1",
                 FE = c("iso3c", "year"),
                estimator = "binning",
                wald = F,
                vcov.type = "cluster",
                cl = "iso3c",
                data = analysis_df, na.rm=T)
ga11 <- interflex(Y = "scaled_cf_principal_received",
                  D = "annual_average_temp_lag1", 
                  X = "gdp_pcmkt", 
                  Z = "annual_W_lag1",
                  FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
ga10 <- interflex(Y = "scaled_cf_total_received",
                  D = "annual_average_temp_lag1", 
                  X = "gdp_pcmkt", 
                  Z = "annual_W_lag1",
                  FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)

## Hottest season temperature * GDP per capita (logged)
gh1 <- interflex(Y = "scaled_climate_laws",
                 D = "hottest_season_temp_lag1", 
                 X = "gdp_pcmkt", 
                 Z = "season_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
gh2 <- interflex(Y = "scaled_ecp_zeros",
                 D = "hottest_season_temp_lag1", 
                 X = "gdp_pcmkt", 
                 Z = "season_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
gh3 <- interflex(Y = "scaled_elec_renew",
                 D = "hottest_season_temp_lag1", 
                 X = "gdp_pcmkt", 
                 Z = "season_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
gh4 <- interflex(Y = "fit",
                 D = "hottest_season_temp_lag1", 
                 X = "gdp_pcmkt", 
                 Z = "season_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
gh5 <- interflex(Y = "quota",
                 D = "hottest_season_temp_lag1", 
                 X = "gdp_pcmkt", 
                 Z = "season_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
gh6 <- interflex(Y = "scaled_delegates",
                 D = "hottest_season_temp_lag1", 
                 X = "gdp_pcmkt", 
                 Z = "season_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
gh7 <- interflex(Y = "scaled_gcg",
                 D = "hottest_season_temp_lag1", 
                 X = "gdp_pcmkt", 
                 Z = "season_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
gh9 <- interflex(Y = "scaled_cf_principal_provided",
                 D = "hottest_season_temp_lag1", 
                 X = "gdp_pcmkt", 
                 Z = "season_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
gh8 <- interflex(Y = "scaled_cf_total_provided",
                 D = "hottest_season_temp_lag1", 
                 X = "gdp_pcmkt", 
                 Z = "season_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
gh11 <- interflex(Y = "scaled_cf_principal_received",
                  D = "hottest_season_temp_lag1", 
                  X = "gdp_pcmkt", 
                  Z = "season_W_lag1",
                  FE = c("iso3c", "year"),
                  estimator = "binning",
                  wald = F,
                  vcov.type = "cluster",
                  cl = "iso3c",
                  data = analysis_df, na.rm=T)
gh10 <- interflex(Y = "scaled_cf_total_received",
                  D = "hottest_season_temp_lag1", 
                  X = "gdp_pcmkt", 
                  Z = "season_W_lag1",
                  FE = c("iso3c", "year"),
                  estimator = "binning",
                  wald = F,
                  vcov.type = "cluster",
                  cl = "iso3c",
                  data = analysis_df, na.rm=T)

## Annual average temperature * VDEM democracy
da1 <- interflex(Y = "scaled_climate_laws",
                 D = "annual_average_temp_lag1", 
                 X = "v2x_polyarchy", 
                 Z = "annual_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
da2 <- interflex(Y = "scaled_ecp_zeros",
                 D = "annual_average_temp_lag1", 
                 X = "v2x_polyarchy", 
                 Z = "annual_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
da3 <- interflex(Y = "scaled_elec_renew",
                 D = "annual_average_temp_lag1", 
                 X = "v2x_polyarchy", 
                 Z = "annual_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
da4 <- interflex(Y = "fit",
                 D = "annual_average_temp_lag1", 
                 X = "v2x_polyarchy", 
                 Z = "annual_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
da5 <- interflex(Y = "quota",
                 D = "annual_average_temp_lag1", 
                 X = "v2x_polyarchy", 
                 Z = "annual_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
da6 <- interflex(Y = "scaled_delegates",
                 D = "annual_average_temp_lag1", 
                 X = "v2x_polyarchy", 
                 Z = "annual_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
da7 <- interflex(Y = "scaled_gcg",
                 D = "annual_average_temp_lag1", 
                 X = "v2x_polyarchy", 
                 Z = "annual_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
da9 <- interflex(Y = "scaled_cf_principal_provided",
                 D = "annual_average_temp_lag1", 
                 X = "v2x_polyarchy", 
                 Z = "annual_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
da8 <- interflex(Y = "scaled_cf_total_provided",
                 D = "annual_average_temp_lag1", 
                 X = "v2x_polyarchy", 
                 Z = "annual_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
da11 <- interflex(Y = "scaled_cf_principal_received",
                  D = "annual_average_temp_lag1", 
                  X = "v2x_polyarchy", 
                  Z = "annual_W_lag1",
                  FE = c("iso3c", "year"),
                  estimator = "binning",
                  wald = F,
                  vcov.type = "cluster",
                  cl = "iso3c",
                  data = analysis_df, na.rm=T)
da10 <- interflex(Y = "scaled_cf_total_received",
                  D = "annual_average_temp_lag1", 
                  X = "v2x_polyarchy", 
                  Z = "annual_W_lag1",
                  FE = c("iso3c", "year"),
                  estimator = "binning",
                  wald = F,
                  vcov.type = "cluster",
                  cl = "iso3c",
                  data = analysis_df, na.rm=T)

## Hottest season temperature * VDEM democracy
dh1 <- interflex(Y = "scaled_climate_laws",
                 D = "hottest_season_temp_lag1", 
                 X = "v2x_polyarchy", 
                 Z = "season_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
dh2 <- interflex(Y = "scaled_ecp_zeros",
                 D = "hottest_season_temp_lag1", 
                 X = "v2x_polyarchy", 
                 Z = "season_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
dh3 <- interflex(Y = "scaled_elec_renew",
                 D = "hottest_season_temp_lag1", 
                 X = "v2x_polyarchy", 
                 Z = "season_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
dh4 <- interflex(Y = "fit",
                 D = "hottest_season_temp_lag1", 
                 X = "v2x_polyarchy", 
                 Z = "season_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
dh5 <- interflex(Y = "quota",
                 D = "hottest_season_temp_lag1", 
                 X = "v2x_polyarchy", 
                 Z = "season_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
dh6 <- interflex(Y = "scaled_delegates",
                 D = "hottest_season_temp_lag1", 
                 X = "v2x_polyarchy", 
                 Z = "season_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
dh7 <- interflex(Y = "scaled_gcg",
                 D = "hottest_season_temp_lag1", 
                 X = "v2x_polyarchy", 
                 Z = "season_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
dh9 <- interflex(Y = "scaled_cf_principal_provided",
                 D = "hottest_season_temp_lag1", 
                 X = "v2x_polyarchy", 
                 Z = "season_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
dh8 <- interflex(Y = "scaled_cf_total_provided",
                 D = "hottest_season_temp_lag1", 
                 X = "v2x_polyarchy", 
                 Z = "season_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
dh11 <- interflex(Y = "scaled_cf_principal_received",
                  D = "hottest_season_temp_lag1", 
                  X = "v2x_polyarchy", 
                  Z = "season_W_lag1",
                  FE = c("iso3c", "year"),
                  estimator = "binning",
                  wald = F,
                  vcov.type = "cluster",
                  cl = "iso3c",
                  data = analysis_df, na.rm=T)
dh10 <- interflex(Y = "scaled_cf_total_received",
                  D = "hottest_season_temp_lag1", 
                  X = "v2x_polyarchy", 
                  Z = "season_W_lag1",
                  FE = c("iso3c", "year"),
                  estimator = "binning",
                  wald = F,
                  vcov.type = "cluster",
                  cl = "iso3c",
                  data = analysis_df, na.rm=T)


gdp_annual <- rbind(ga1$est.bin[[1]], ga2$est.bin[[1]],
      ga3$est.bin[[1]], ga4$est.bin[[1]],
      ga5$est.bin[[1]], ga6$est.bin[[1]],
      ga7$est.bin[[1]], ga8$est.bin[[1]],
      ga9$est.bin[[1]], ga10$est.bin[[1]], 
      ga11$est.bin[[1]]) %>% 
  as.data.frame() %>% 
  select(beta = coef, se = sd, ci_lower = CI.lower, ci_upper = CI.upper) %>% 
  mutate(ci_lower_correction = beta - 2.84*se,
         ci_upper_correction = beta + 2.84*se) %>% 
  mutate(Bin = rep(c("Low", "Medium", "High"), times = 11)) %>% 
  mutate(Bin = factor(Bin, c("Low", "Medium", "High"))) %>% 
  mutate(outcome = rep(select(analysis_df, starts_with("scaled")) %>% names,
                       each = 3)) %>% 
  mutate(moderator = "GDP per capita",
         treatment = "Annual average")

gdp_hottest <- rbind(gh1$est.bin[[1]], gh2$est.bin[[1]],
                    gh3$est.bin[[1]], gh4$est.bin[[1]],
                    gh5$est.bin[[1]], gh6$est.bin[[1]],
                    gh7$est.bin[[1]], gh8$est.bin[[1]],
                    gh9$est.bin[[1]], gh10$est.bin[[1]], 
                    gh11$est.bin[[1]]) %>% 
  as.data.frame() %>% 
  select(beta = coef, se = sd, ci_lower = CI.lower, ci_upper = CI.upper) %>% 
  mutate(ci_lower_correction = beta - 2.84*se,
         ci_upper_correction = beta + 2.84*se) %>% 
  mutate(Bin = rep(c("Low", "Medium", "High"), times = 11)) %>% 
  mutate(Bin = factor(Bin, c("Low", "Medium", "High"))) %>% 
  mutate(outcome = rep(select(analysis_df, starts_with("scaled")) %>% names,
                       each = 3)) %>% 
  mutate(moderator = "GDP per capita",
         treatment = "Hottest season")

dem_annual <- rbind(da1$est.bin[[1]], da2$est.bin[[1]],
                    da3$est.bin[[1]], da4$est.bin[[1]],
                    da5$est.bin[[1]], da6$est.bin[[1]],
                    da7$est.bin[[1]], da8$est.bin[[1]],
                    da9$est.bin[[1]], da10$est.bin[[1]], 
                    da11$est.bin[[1]]) %>% 
  as.data.frame() %>% 
  select(beta = coef, se = sd, ci_lower = CI.lower, ci_upper = CI.upper) %>% 
  mutate(ci_lower_correction = beta - 2.84*se,
         ci_upper_correction = beta + 2.84*se) %>% 
  mutate(Bin = rep(c("Low", "Medium", "High"), times = 11)) %>% 
  mutate(Bin = factor(Bin, c("Low", "Medium", "High"))) %>% 
  mutate(outcome = rep(select(analysis_df, starts_with("scaled")) %>% names,
                       each = 3)) %>% 
  mutate(moderator = "Democracy",
         treatment = "Annual average")

dem_hottest <- rbind(dh1$est.bin[[1]], dh2$est.bin[[1]],
                     dh3$est.bin[[1]], dh4$est.bin[[1]],
                     dh5$est.bin[[1]], dh6$est.bin[[1]],
                     dh7$est.bin[[1]], dh8$est.bin[[1]],
                     dh9$est.bin[[1]], dh10$est.bin[[1]], 
                     dh11$est.bin[[1]]) %>% 
  as.data.frame() %>% 
  select(beta = coef, se = sd, ci_lower = CI.lower, ci_upper = CI.upper) %>% 
  mutate(ci_lower_correction = beta - 2.84*se,
         ci_upper_correction = beta + 2.84*se) %>% 
  mutate(Bin = rep(c("Low", "Medium", "High"), times = 11)) %>% 
  mutate(Bin = factor(Bin, c("Low", "Medium", "High"))) %>% 
  mutate(outcome = rep(select(analysis_df, starts_with("scaled")) %>% names,
                       each = 3)) %>% 
  mutate(moderator = "Democracy",
         treatment = "Hottest season")


# Plot in four facets
rbind.data.frame(gdp_annual, gdp_hottest, 
                 dem_annual, dem_hottest) %>% 
  remove_rownames() %>% 
  group_by(Bin, moderator, treatment) %>% 
  mutate(n = row_number()) %>% 
  group_by(outcome, moderator, treatment) %>% 
  mutate(m = row_number()) %>% 
  ungroup() %>% 
  mutate(index = ifelse(m == 1, n-0.2,
                        ifelse(m==3, n+0.2, n))) %>% 
  # Create facet ID
  mutate(Facet = paste(treatment, moderator, sep = " * ")) %>% 
  ggplot(., aes(beta, index)) +
  facet_wrap(~ Facet, nrow = 2) +
  coord_cartesian(xlim = c(-0.45,0.45)) +
  scale_x_continuous(breaks = c(-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4)) +
  geom_segment(aes(x = ci_lower_correction, xend = ci_upper_correction,
                   y = index, yend = index),
               size = 1, color = "#79b321") +
  geom_segment(aes(x = ci_lower, xend = ci_upper, 
                   y = index, yend = index), 
               size = 2, color = "#456613") +
  geom_point(size = 3.5, aes(shape = Bin)) +
  scale_shape_manual(breaks=c("Low","Medium","High"),
                     values = c("L","M","H")) + # 16, 3, 4
  scale_color_manual(values = c("black", "black", "black")) +
  geom_vline(xintercept = 0, linetype = 1, colour = "black") +
  scale_y_continuous(trans = "reverse",
                     breaks = c(1:11),
                     labels = unique(labels_vector$labels)) +
  labs(x="Regression coefficient and confidence intervals", 
       y="") +
  geom_vline(xintercept = 0) +
  theme(panel.spacing = unit(1, "lines")) +
  theme_tufte(base_family = "Helvetica") +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill=NA),
        panel.grid.major.x = element_line(colour = "grey90"),
        panel.grid.minor.x = element_line(colour = "grey90"),
        # panel.grid.major.y = element_line(colour = "grey90"),
        panel.ontop = FALSE,
        text=element_text(size=16))
ggsave("text/coefplot_heterogeneous_continuous.pdf", units = "in", height= 10, width = 10)

  
